
# Set ylims on divergence, number of clones
ylimit <- NULL
ylimit_clones <- NULL
ylimit[1:3] <- 0.15
ylimit_clones[1:3] <- 2000
ylimit[4:6] <- 0.05
ylimit_clones[4:6] <- 2000


for(i in 1:1){
	png(paste("sweep",i,"/sweep",i,".divergence.png",sep=""),width=4.1,height=4.1,units="in",res=72)
	par(mfrow=c(1,1))
	par(pty="s")
	par(mar=c(3,3,3,3),cex=0.8)
	divergence <- read.csv(paste("sweep",i,"/sweep",i,".divergence.log",sep=""),header=FALSE)
	plot(divergence[,1]/364,divergence[,2],xlab=NA,ylab=NA,t="l",ylim=c(0,ylimit[i]),xlim=c(0,7300/364),col="black")
	fixations <- read.csv(paste("sweep",i,"/sweep",i,".fixations.log",sep=""),header=FALSE)
	fixations <- cbind(fixations[,1:4], divergence[,2])
	fixev <- fixations[fixations[,3]>0,]	
	points(fixev[,1]/364,fixev[,5]-ylimit[i]/30,col="black",pch=fixev[,3]+14,cex=1.5)
	fixev <- fixations[fixations[,4]>0,]
	points(fixev[,1]/364,fixev[,5]+ylimit[i]/30,col="red",pch=fixev[,4]+14,cex=1.5)
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1]/364,fixev[,5],col="black",pch="|",cex=1.0)

	# Add on number of clones
	clones_total <- read.csv(paste("sweep",i,"/sweep",i,".unique.clones.log",sep=""),header=FALSE)
	# Re-scale the clones number data between 0.00 and 0.04
	rescaled <- (clones_total[,2]/ylimit_clones[i])*ylimit[i]
	lines(clones_total[,1]/364, rescaled, col="black",lty=2, ylab="Number of clones")
	axis(side=4,at=ylimit[i]*c(0:5)/5,labels=ylimit_clones[i]*c(0:5)/5)
	f <- divergence[nrow(divergence)/4*c(1:5),3]
	axis(side=3,at=20*c(1:5)/4,labels=f)

	mtext("Years", side=1, cex=par("cex"), line = 2)
	mtext("Divergence", side=2, cex=par("cex"), line = 2)
	mtext("Generations", side=3, cex=par("cex"), line = 2)
	mtext("Number of clones", side=4, cex=par("cex"), line = 2)
	dev.off()
}

# Plot frequency over time for all clones that have reached more than 5% frequency any time throughout 20 years
for(i in 1:1){
	pdf(paste("sweep",i,"/sweep",i,".frequency.pdf",sep=""),width=2,height=2,pointsize=8)
	par(mfrow=c(1,1))
	par(pty="s")
	par(mar=c(3,3,1,1))
	par(mgp=c(2,1,0))
	d <- read.csv(paste("sweep",i,"/sweep",i,".frequencies",sep=""),header=FALSE)
	dat <- d[d[,3]>0.05&d[,1]>0,]
	if(nrow(dat)!=0){
	# get the frequency of clones in dat, even if it dropped to 0 at some time points
	clones.freq <- d[d[,1] %in% dat[,1],]
	plot(clones.freq[,2]/364,clones.freq[,3],col="black",ylab=NA,xlab=NA,pch=NA,ylim=c(0,1),xlim=c(0,20))
	rrr <- unique(clones.freq[,1])
	for(i in 1:length(rrr)){
		lines(clones.freq[clones.freq[,1]==rrr[i],2]/364,clones.freq[clones.freq[,1]==rrr[i],3],col="black")
	}
	} else {
	plot(0,0,col="black",ylab=NA,xlab=NA,pch=NA,ylim=c(0,1),xlim=c(0,20))
	}
	mtext("Years", side=1, cex=par("cex"), line = 2)
	mtext("Frequency", side=2, cex=par("cex"), line = 2)
	dev.off()
}






# Plot clone size distribution over time
par(mfrow=c(1,1))
for(i in 2:2){
	sizedistr <- read.csv(paste("sweep",i,"/sweep",i,".unique.clones.size.distrib.log",sep=""))
	boxplot(split(sizedistr[,2], sizedistr[,1]),ylim=c(0,200))
	clonenums <- as.numeric(sapply(split(sizedistr[,2], sizedistr[,1]),length))
}








# Set ylim on divergence
ylimit <- 0.04
# s=0.5, selmu = 10^-6, 10^-7, 10^-8 # s=1, selmu = 10^-6, 10^-7, 10^-8
par(mfrow=c(2,3))
for(i in 7:8){
	divergence <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".divergence.log",sep=""))
	plot(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="black")
	fixations <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".fixations.log",sep=""))
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="black",pch=fixev[,5]+14,cex=1.5)

	divergence <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="red")
	fixations <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".fixations.log",sep=""))
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="red",pch=fixev[,5]+14,cex=1.5)

	divergence <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="blue")
	fixations <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".fixations.log",sep=""))
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="blue",pch=fixev[,5]+14,cex=1.5)
}

# Repro driven
# Set ylim on divergence
ylimit <- 0.04
# s=0.5, selmu = 10^-6, 10^-7, 10^-8 # s=1, selmu = 10^-6, 10^-7, 10^-8
par(mfrow=c(2,3))
for(i in 9:12){
	divergence <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".divergence.log",sep=""))
	plot(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="black")
	fixations <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".fixations.log",sep=""))
	divergence <- divergence[1:nrow(fixations),]
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="black",pch=fixev[,5]+14,cex=1.5)

	divergence <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="red")
	fixations <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".fixations.log",sep=""))
	divergence <- divergence[1:nrow(fixations),]
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="red",pch=fixev[,5]+14,cex=1.5)

	divergence <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="blue")
	fixations <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".fixations.log",sep=""))
	divergence <- divergence[1:nrow(fixations),]
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="blue",pch=fixev[,5]+14,cex=1.5)
}

# Repro driven
# Set ylim on divergence
ylimit <- 0.04
# s=0.5, selmu = 10^-6, 10^-7, 10^-8 # s=1, selmu = 10^-6, 10^-7, 10^-8
#par(mfrow=c(2,3))
for(i in 15:15){
	divergence <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".divergence.log",sep=""))
	plot(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="black")
	fixations <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".fixations.log",sep=""))
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="black",pch=fixev[,5]+14,cex=1.5)

	divergence <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="red")
	fixations <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".fixations.log",sep=""))
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="red",pch=fixev[,5]+14,cex=1.5)

	divergence <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="blue")
	fixations <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".fixations.log",sep=""))
	fixations <- cbind(fixations, divergence[,2])
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],fixev[,3],col="blue",pch=fixev[,5]+14,cex=1.5)
}









for(i in 3:5){
	divergence <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".divergence.log",sep=""))
	plot(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="black")
	fixations <- read.csv(paste("sweep",i*3+1,"/sweep",i*3+1,".fixations.log",sep=""))
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],divergence[,2],col=fixev[,2])

	divergence <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="red")
	fixations <- read.csv(paste("sweep",i*3+2,"/sweep",i*3+2,".fixations.log",sep=""))
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],divergence[,2],col=fixev[,2])

	divergence <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".divergence.log",sep=""))
	lines(divergence[,1],divergence[,2],xlab="Time",ylab="Divergence",t="l",ylim=c(0,ylimit),xlim=c(182,7300),col="blue")
	fixations <- read.csv(paste("sweep",i*3+3,"/sweep",i*3+3,".fixations.log",sep=""))
	fixev <- fixations[fixations[,2]>0,]
	points(fixev[,1],divergence[,2],col=fixev[,2])
}




d <- read.csv("sweep21/sweep21.frequencies",header=FALSE)

dat <- d[d[,3]>0.05&d[,1]>0,]

rrr <- unique(dat[,1])
rrrc <- gray(rrr/max(rrr))
clonecol <- cbind(rrr,rrrc)

plot(dat[,2],dat[,3],col="black",ylab="Frequency",xlab="Time (days)",pch=19)
for(i in 1:length(rrr)){
lines(dat[dat[,1]==rrr[i],2],dat[dat[,1]==rrr[i],3],col="black")
}

# Color by clone selective advantage

# Load clones selective advantage, color clones by advantage
clones.fitness <- read.csv("sweep21/sweep21.clones",header=FALSE)
d <- read.csv("sweep21/sweep21.frequencies",header=FALSE)

# Remove duplicates
clones.fitness <- unique(clones.fitness)

cfit <- cbind(rep(0,nrow(clones.fitness)),rep(0,nrow(clones.fitness)))
cfit[clones.fitness[,2]==1,1] <- "black"
cfit[clones.fitness[,2]==2,1] <- "orange"
cfit[clones.fitness[,2]==4,1] <- "red"
cfit[clones.fitness[,2]==8,1] <- "brown"
cfit[clones.fitness[,3]==1,2] <- "black"
cfit[clones.fitness[,3]==2,2] <- "orange"
cfit[clones.fitness[,3]==4,2] <- "red"
cfit[clones.fitness[,3]==8,2] <- "brown"

# get clones that have reached more than 5% frequency over 20 years
# clone_id, time, frequency
dat <- d[d[,3]>0.05&d[,1]>0,]

intersect <- function(x, y) y[match(x, y, nomatch = 0)]

# get the frequency of clones in dat, even if it dropped to 0 at some time points
clones.freq <- d[d[,1] %in% dat[,1],]

# Plot
par(mfrow=c(1,2))
plot(clones.freq[,2],clones.freq[,3],col="black",ylab="Frequency",xlab="Time (days)",main="Survival Advantage",pch=NA)
rrr <- unique(clones.freq[,1])
for(i in 1:length(rrr)){
lines(clones.freq[clones.freq[,1]==rrr[i],2],clones.freq[clones.freq[,1]==rrr[i],3],col=cfit[clones.fitness[,1]==rrr[i],1])
}
plot(clones.freq[,2],clones.freq[,3],col="black",ylab="Frequency",xlab="Time (days)",main="Reproductive Advantage",pch=NA)
rrr <- unique(clones.freq[,1])
for(i in 1:length(rrr)){
lines(clones.freq[clones.freq[,1]==rrr[i],2],clones.freq[clones.freq[,1]==rrr[i],3],col=cfit[clones.fitness[,1]==rrr[i],2])
}



# Load data
clonesn21 <- read.csv("sweep21/sweep21.clones.number.log",header=FALSE)

par(mfrow=c(1,1))
plot(clonesn21[,1],clonesn21[,2],col="black",ylab="Absolute number of clones",xlab="Time (days)",main="10e-06",type="l")




# Plot number of clones over time
# Load data
clonesn6 <- read.csv("sweep3/sweep3.clones.number.log",header=FALSE)
clonesn7 <- read.csv("sweep5/sweep5.clones.number.log",header=FALSE)
clonesn8 <- read.csv("sweep6/sweep6.clones.number.log",header=FALSE)

par(mfrow=c(1,3))
plot(clonesn6[,1],clonesn6[,2],col="black",ylab="Absolute number of clones",xlab="Time (days)",main="10e-06",type="l")
plot(clonesn7[,1],clonesn7[,2],col="black",ylab="Absolute number of clones",xlab="Time (days)",main="10e-07",type="l")
plot(clonesn8[,1],clonesn8[,2],col="black",ylab="Absolute number of clones",xlab="Time (days)",main="10e-08",type="l")

grid6 <- read.csv("sweep6/sweep6.log",header=FALSE)
ggrid6 <- grid6[,2:3]
cnum <- NULL
ctime <- NULL
for(i in 1:40){
	cnum[i] <- nrow(unique(ggrid6[ggrid6[,2]==i*182,1:2]))
	ctime[i] <- i*182
}
par(mfrow=c(1,1))
plot(ctime,cnum,col="black",ylab="Absolute number of clones in grid",xlab="Time (days)",main="10e-08",type="l")







# Average cell lifespan
distrib <- NULL
for(i in 1:1000){

	lifespan <- 0	
	div <- rexp(1,1/4)
	death <- rexp(1,1/40)
	while(div<death){
		lifespan <- lifespan + 1
		div <- rexp(1,1/4)
		death <- rexp(1,1/40)
	}	
	distrib[i] <- lifespan
}
hist(distrib,freq=FALSE)






